* Required package - rdrobust
*net install st0366_1.pkg

* Restricted relication data
u restricteddata, clear


*********************************************************************************
       **************************** FIGURES *****************************
*********************************************************************************

* graph format
set scheme s1mono
graph set window fontface "Times"



*******************************************
* Figure 1
*******************************************
* (select and run everything through "graph combine")


* Creating granular bins for RD graphs
cap g k=.
qui forv j=-3(.05)3 {
	replace k=`j' if inrange( std_admscore ,`j',`j'+.05) 
}
bys k: egen gbmd = mean( std_admscore ) 
bys k : egen gad = mean(enrolled) 
bys k : g n = _n==1 

* 1st-Stage RD estimate
rdrobust enrolled std_admscore, vce(cluster id) all 
loc tau = string(e(tau_bc),"%9.3f")
loc star = cond(e(pv_rb)<=.01,"***",cond(e(pv_rb)<=.05,"**", ///
	cond(e(pv_rb)<.1,"*","")))
loc se = string(e(se_tau_rb),"%9.3f")

tw sc gad gbmd if std_admsc>-1 & std_admsc<1 & n==1, mcolor(gs9) ///
		msize(medium) msymbol(circle) || ///
	qfitci enrolled std_admscore if std_admsc>-1 & std_admsc<0, ///
		lcolor(blue) alcolor(blue%8) acolor( blue%8) || ///
	qfitci enrolled std_admscore if std_admsc>=0 & std_admsc<1, ///
		lcolor(blue) alcolor(blue%8) acolor( blue%8) || ///
	, xline(0, lpattern(dash) lcolor(red)) leg(off) ///
		text(.3 0 "{&tau} = `tau'`star'" "   (`se')", place(e) size(medium)) ///
		xlab(0 -1(.2)1 ) xti(Admission score) ///
		ylabel(, nogrid) yti(Probability of enrollment) ///
		ti("All sample", color(black)) ///
		graphregion(fcolor(white) color(white)) ///
		yla(0(.2)1) name(enrollment, replace) nodraw


* RD estimate per income group

qui foreach y of varlist high80 low20 {

	* Granular bins
	
	if "`y'"=="high80" loc h .05 
    if "`y'"=="low20" loc h .05 
	
	g k`y'=.
	forv j=-3(`h')3 {
		replace k`y'=`j' if inrange( std_admscore ,`j',`j'+`h') 
	}
	bys k`y': egen gbmd`y' = mean( std_admscore ) if `y'==1
	so k`y' gbmd`y' 
	bys k`y': replace gbmd`y' = gbmd`y'[1] if `y'==1
	bys k`y': egen gad`y' = mean(enrolled) if `y'==1
	so k`y' gad`y' 
	bys k`y': replace gad`y' = gad`y'[1] if `y'==1
	bys k`y': g n`y' = _n==1 if `y'==1

	
	* Graph titles
	if "`y'"=="high80" loc mytitle "High income"
    if "`y'"=="low20" loc mytitle "Low income"

	
	* 1st-stage RD estimate
	rdrobust enrolled std_admscore if `y'==1, vce(cluster id) all 

    loc tau = string(e(tau_bc),"%9.3f")
    loc star = cond(e(pv_rb)<=.01,"***",cond(e(pv_rb)<=.05,"**", ///
		cond(e(pv_rb)<.1,"*","")))
    loc se = string(e(se_tau_rb),"%9.3f")

	tw sc gad`y' gbmd`y' if std_admsc>-1 & std_admsc<1 & n`y'==1 & `y'==1, ///
			mcolor(gs9) msize(medium) msymbol(circle) || ///
		qfitci enrolled std_admscore if std_admsc>-1 & std_admsc<=0 & `y'==1, ///
			lcolor(blue) alcolor(blue%8) acolor( blue%8) || ///
		qfitci enrolled std_admscore if std_admsc>0 & std_admsc<1 & `y'==1, ///
			lcolor(blue) alcolor(blue%8) acolor( blue%8) || ///
		, xline(0, lpattern(dash) lcolor(red)) leg(off) ///
			text(.3 0 "{&tau} = `tau'`star'" "   (`se')", place(e) size(medium)) ///
			xlab(0 -1(.2)1 ) xti(Admission score) ylabel(, nogrid) ///
			ti("`mytitle'", color(black)) graphregion(fcolor(white) color(white)) ///
			yla(0(.2)1) yti(Probability of enrollment) ///
			name(enrollment`y', replace) nodraw
		drop k`y' gbmd`y' gad`y' n`y'
}

* Final figure
graph combine enrollment enrollmentlow20 enrollmenthigh80, r(1) ///
	graphregion(color(white)) xsize(12) iscale(1) ///
    name(figure_income, replace) 

	

*******************************************
* Figure 2 
*******************************************
* (select and run everything through "graph combine")


* Outcome
loc r "crime_adm10"

* Generate graph per income group

qui foreach y in all high80 low20 {

	* Graph titles
	if "`y'"=="high80" loc f "High income"
    if "`y'"=="low20" loc f "Low income"
    if "`y'"=="all" loc f "All sample"

	* Graph parameters
    loc w "-.03(.03).12"
	loc z 0.09 


	* Granular bins
	g k`y'=.
	forv j=-3(.05)3 {
		replace k`y'=`j' if inrange( std_admscore ,`j',`j'+.05) 
	}

	bys k`y': egen gbmd`y' = mean( std_admscore ) if `y'==1 
	so k`y' gbmd`y' 
	bys k`y': replace gbmd`y' = gbmd`y'[1] if `y'==1
	bys k`y': egen gad`y' = mean(`r') if `y'==1 
	so k`y' gad`y' 
	bys k`y': replace gad`y' = gad`y'[1] if `y'==1
	bys k`y': g n`y' = _n==1 
	
	
	* RD estimate
	rdrobust `r' std_admscore if `y'==1 , vce(cluster id) all 

    loc tau = string(e(tau_bc),"%9.3f")
    loc star = cond(e(pv_rb)<=.01,"***",cond(e(pv_rb)<=.05,"**", ///
		cond(e(pv_rb)<.1,"*","")))
    loc se = string(e(se_tau_rb),"%9.3f")

	tw sc gad`y' gbmd`y' if std_admsc>-1 & std_admsc<1 & n`y'==1 & `y'==1, ///
			mcolor(gs9) msize(medium) msymbol(circle) || ///
		lfitci `r' std_admscore if std_admsc>-1 & std_admsc<=0 & `y'==1 , ///
			lcolor(blue) alcolor(blue%8) acolor( blue%8) level(95) || ///
		lfitci `r' std_admscore if std_admsc>=0 & std_admsc<1 & `y'==1 , ///
			lcolor(blue) alcolor(blue%8) acolor( blue%8) level(95) || ///
		, xline(0, lpattern(dash) lcolor(red)) leg(off) ///
			text("`z'" 0 "{&tau} = `tau'`star'" "  (`se')", place(e) size(medium)) ///
			xlab(0 -1(.2)1 ) xti(Admission score) ylabel(, nogrid) ///
			ti("`f'", color(black)) ylabel(`w', nogrid) ///
			graphregion(fcolor(white) color(white)) name(`r'`y', replace) nodraw
	drop k`y' gbmd`y' gad`y' n`y' 
}		

* Final figure
graph combine `r'all `r'low20 `r'high80, r(1) graphregion(color(white)) ///
	xsize(12) iscale(1) name(`r', replace)	



********************************************
* Figure 3 
********************************************

**********************
* Figure 3 (Panel A)
**********************
* (select and run everything through "graph combine")

* Types of crime

loc ys
foreach y in crime_adm10 crime_violent_adm10 crime_property_adm10 ///
	crime_drugs_adm10 crime_public_adm10 crime_public_drug_adm10 ///
	crime_traffic_adm10 crime_unclassified_adm10 {

	loc ys "`ys' `y'"
}


* Graphs per income group
qui foreach g in all low20 high80 {

	loc gg
	loc ti "All sample"
	loc yl "-.1(.03).02"
	if "`g'"=="low20" {
		loc gg "l"
		loc ti "Low income"
		loc yl "-.16(.04).04"
	}
	if "`g'"=="high80" {
		loc gg "h"
		loc ti "High income"
	}

	* Location of estimates in the horizontal axis
	g double `gg'r_m=. // males
	g double `gg'r_f=. // females

	* Variables to store estimates
		* Males
	g `gg'up90_m=.
	g `gg'low90_m=.
    g `gg'up95_m=.
	g `gg'low95_m=.
	g `gg'coef_m=.

		* Females
	g `gg'up90_f=.
	g `gg'low90_f=.
    g `gg'up95_f=.
	g `gg'low95_f=.
	g `gg'coef_f=.		


	* Effect per type of crime

	local z=1
	foreach y in `ys' {

		replace `gg'r_m = round(`z'-.1, .0) in `z'
		replace `gg'r_f = round(`z'+.1, .0) in `z'

	
		* Effect per gender

		foreach s in male female {
		
			if "`s'"=="male" 	loc p "m"
			if "`s'"=="female" 	loc p "f"
			
			* RD estimate

			cap rdrobust `y' std_admscore if `s'==1 & `g'==1, ///
				vce(cluster id) all
		
			if _rc==0 {
				replace `gg'up90_`p' = e(tau_bc) + 1.645*e(se_tau_rb) in `z'
				replace `gg'low90_`p' = e(tau_bc) - 1.645*e(se_tau_rb) in `z'
				replace `gg'up95_`p' = e(tau_bc) + 1.96*e(se_tau_rb) in `z'
				replace `gg'low95_`p' = e(tau_bc) - 1.96*e(se_tau_rb) in `z'
				replace `gg'coef_`p' = e(tau_bc) in `z'
			}
			else {
				su `y' if `s'==1 & `g'==1 & inrange(std_admscore,-2,2)
				loc rm = r(mean)
				if `rm' < 10e-4 {
					replace `gg'up90_`p' = 0 in `z'
					replace `gg'low90_`p' = 0 in `z'
					replace `gg'up95_`p' = 0 in `z'
					replace `gg'low95_`p' = 0 in `z'
					replace `gg'coef_`p' = 0 in `z'
				}
				else {
					di "Error, `y', `s', `g'"
				}
			}
		}
		loc z = `z'+1
	}
	twoway (rcap `gg'up95_m `gg'low95_m `gg'r_m, lcolor(blue) ///
				msize(*0.7) color(blue)) ///
	       (rcap `gg'up90_m `gg'low90_m `gg'r_m, lcolor(blue) ///
				msize(*1.3) color(blue)) ///
		   (scatter `gg'coef_m `gg'r_m, mcolor(blue) msymbol(circle) ///
				msize(medium)) ///
		   (rcap `gg'up95_f `gg'low95_f `gg'r_f, lcolor(gray) ///
				msize(*0.7) color(gray)) ///
	       (rcap `gg'up90_f `gg'low90_f `gg'r_f, lcolor(gray) ///
				msize(*1.3) color(gray)) ///
		   (scatter `gg'coef_f `gg'r_f, mcolor(gray) msymbol(diamond) ///
				msize(medium)) ///	
		, yline(0, lpattern(dash) lcolor(red)) ylabel(`yl', nogrid) ///
            graphregion(fcolor(white) color(white)) xti("") /// 
			title(`ti', color(black) ) xla(1 `" "All" "crimes" "' ///
				2 `" "Violent" "crimes" "' 3 `" "Property" "crimes" "' ///
				4 `" "Drug-" "related" "crimes" "' ///
				5 `" "Against" "public" "interest" "' ///
				6 `" "Against" "public" "+ drugs" "' ///
				7 `" "Traffic-" "related" "crimes" "' ///
				8 `" "Unclassif." "crimes" "', angle(75) labsize(small)) ///
			xsize(8) yti(Estimated discontinuity) ///
			legend(pos(5) ring(0) col(2) rows(1) size(3) ///
			order(3 "Male" 6 "Female" ) region(style(none)) ) ///
			name(`gg'gender, replace) nodraw
}


* Final figure (panel A)
gr combine gender lgender hgender, r(1) graphregion(color(white)) xsize(15) ///
	iscale(1.3) name(gendert, replace)


drop r_m-hcoef_f


**********************	
* Figure 3 (Panel B)
**********************	
* (select and run everything through "graph combine")

* Types of crime

loc ys
foreach y in crime_adm10 crime_violent_adm10 crime_property_adm10 ///
	crime_drugs_adm10 crime_public_adm10 crime_public_drug_adm10 ///
	crime_traffic_adm10 crime_unclassified_adm10 {

	loc ys "`ys' `y'"
}


* Bandwidth variable
cap g kw = .


* Graphs per income group
qui foreach g in all low20 high80 {

	loc gg
	loc ti "All sample"
	if "`g'"=="low20" {
		loc gg "l"
		loc ti "Low income"
	}
	if "`g'"=="high80" {
		loc gg "h"
		loc ti "High income"
	}

	* Location of estimates in the horizontal axis
	g double `gg'r_m=. // males
	g double `gg'r_f=. // females

	* Variables to store estimates
	g `gg'coef_m=. // males
	g `gg'coef_f=. // females

	
	* Estimate per type of crime
	
	local z=1
	foreach y in `ys' {
		
		replace `gg'r_m = round(`z'-.1,.0) in `z'
		replace `gg'r_f = round(`z'+.1,.0) in `z'


		* Estimate per gender

		foreach s in male female {
		
			if "`s'"=="male" 	loc p "m"
			if "`s'"=="female" 	loc p "f"

			* Optimal bandwidth			
			rdbwselect `y' std_admscore if `s'==1 & `g'==1, vce(cluster id)
			loc bw = e(h_mserd)
			if `bw'==. {
				loc bw = 2
			}
			replace kw = max(0,1 - abs(std_admscore)/`bw')


			* Estimation of Y(0) 
			reg `y' if std_admscore<=0 & `s'==1 & `g'==1 [aw=kw] 
			replace `gg'coef_`p' = round(_b[_cons],.001) in `z'
		}
	local z = `z' + 1
	}
	
	twoway (bar `gg'coef_m `gg'r_m, color(blue%20) barwidth(.2)) ///
		   (bar `gg'coef_f `gg'r_f, color(gs10) barwidth(.2)) ///	
		, yline(0, lpattern(dash) lcolor(red)) ylabel(0(.02).1, nogrid) ///
		graphregion(fcolor(white) color(white)) xti("") /// 
		title(`ti', color(black) ) xla(1 `" "All" "crimes" "' ///
			2 `" "Violent" "crimes" "' 3 `" "Property" "crimes" "' ///
			4 `" "Drug-" "related" "crimes" "' 5 `" "Against" "public" "interest" "' ///
			6 `" "Against" "public" "+ drugs" "' 7 `" "Traffic-" "related" "crimes" "' ///
			8 `" "Unclassif." "crimes" "', angle(75) labsize(small)) xsize(8) ///
		yti(Baseline rate) legend(pos(2) ring(0) col(2) rows(1) size(3) ///
		order(1 "Male" 2 "Female" ) region(style(none)) ) ///
		name(`gg'genderb, replace)	nodraw
}

* Final figure (panel B)
gr combine genderb lgenderb hgenderb, r(1) graphregion(color(white)) xsize(15) ///
	iscale(1.3) name(genderbt, replace)


drop r_m-hcoef_f



********************************************
* Figure 4
********************************************
* (select and run everything through "graph combine")

* Bandwidth variable
cap g kw = .

loc xti "Year after application"
loc int "2/12" // time interval
loc w "-.1(.03).02" // y axis label


* Crime per year and cumulative
qui foreach y in crime_adm crime_at {

	* Estimation per income group
	foreach g in all low20 high80 {
		
		loc gg
		loc ti "All sample"
		if "`g'"=="low20" {
			loc gg "l"
			loc ti "Low income"
		}
		if "`g'"=="high80" {
			loc gg "h"
			loc ti "High income"
		}

		* Location of estimates in the horizontal axis
		g r_`y'=.
		
		* Variables storing estimates
		g up90_`y'=.
		g low90_`y'=.
		g up95_`y'=.
		g low95_`y'=.
		g coef_`y'=.	
		g b_`y'=.	

		local z=1
		
		* Estimate per year
		forv i = `int' {

			replace r_`y' = `i' in `z'

			su `y'`i' if `g'==1
			loc mb = r(mean)

			if `mb'>0 {
				
				* Optimal bandwidth
				rdbwselect `y'`i' std_admscore if `g'==1, vce(cluster id)
				loc bw = e(h_mserd)
				replace kw = max(0,1 - abs(std_admscore)/`bw')
			
				* Estimation of baseline rate, Y(0)
				reg `y'`i' std_admscore if std_admscore<=0 & `g'==1 [aw=kw] 
				replace b_`y' = round(_b[_cons],.001) in `z'
				

				* RD estimate
				rdrobust `y'`i' std_admscore if `g'==1, vce(cluster id) all 
		
				replace up90_`y'= e(tau_bc) + 1.645*e(se_tau_rb) in `z'
				replace low90_`y' = e(tau_bc) - 1.645*e(se_tau_rb) in `z'
				replace up95_`y'= e(tau_bc) + 1.96*e(se_tau_rb) in `z'
				replace low95_`y' = e(tau_bc) - 1.96*e(se_tau_rb) in `z'
				replace coef_`y' = e(tau_bc) in `z'			
			}
			else {
				replace up95_`y'=0 in `z'
				replace low95_`y'=0 in `z'
				replace coef_`y'=0 in `z'
				replace up90_`y'=0 in `z'
				replace low90_`y'=0 in `z'
				replace coef_`y'=0 in `z'
				replace b_`y'=0 in `z'
			}
			loc z=`z'+1
		}

		* Graph
		twoway (rarea up95_`y' low95_`y' r_`y', lcolor(ltbluishgray%8) ///
					msize(*0.7) color(blue%8) yaxis(1)) ///
			   (rarea up90_`y' low90_`y' r_`y', lcolor(ltbluishgray%10) ///
					msize(*1.3) color(blue%10) yaxis(1)) ///
			   (connected coef_`y' r_`y', mcolor(gs9) lcolor(blue) ///
					msymbol(circle) msize(medium) yaxis(1)) ///
			   (bar b_`y' r_`y', color(gs8%15) yaxis(2) ///
					ytitle("Baseline rate", axis(2)) ylab(0(.02).08, axis(2))) ///
			   , yline(0, lpattern(dash) lcolor(red)) ylabel(`w', nogrid) ///
					graphregion(fcolor(white) color(white)) /// 
					title(`ti', color(black) ) xla(2/12) xsize(8) ///
					legend(off) yti(Estimated discontinuity) xtitle(`xti') ///
					name(r_`y'`gg', replace) nodraw

		drop up90_* low90_* up95_* low95_* coef_* r_* b_*
	}
}

* Panel A	
graph combine r_crime_at r_crime_atl r_crime_ath, r(1) ///
	graphregion(color(white)) xsize(15) iscale(1.3) name(crime_at, replace)	
	
* Panel B	
graph combine r_crime_adm r_crime_adml r_crime_admh, r(1) ///
	graphregion(color(white)) xsize(15) iscale(1.3) name(crime_adm, replace)



********************************************
* Figure 5
********************************************
* (select and run everything through "graph combine")

* Outcome
loc y "employed_adm"


* Graph parameters
loc xti "Year after application" // x axis title
loc int "0/12" // time interval
loc yl "-.3(.1).2" // y axis label


* Estimates per income group
qui foreach g in all low20 high80 {
		
	loc gg
	loc ti "All sample"
	if "`g'"=="low20" {
		loc gg "l"
		loc ti "Low income"
	}
	if "`g'"=="high80" {
		loc gg "h"
		loc ti "High income"
	}

	* Location of estimates in the horizontal axis
	g r_`y'=.
	
	* Variables storing estimates
	g up90_`y'=.
	g low90_`y'=.
	g up95_`y'=.
	g low95_`y'=.
	g coef_`y'=.	


	* Estimation per year
	local z=1
	forv i = `int' {

		replace r_`y' = `i' in `z'

		* RD estimate
		rdrobust `y'`i' std_admscore if `g'==1, vce(cluster id) all 

		replace up90_`y'= e(tau_bc) + 1.645*e(se_tau_rb) in `z'
		replace low90_`y' = e(tau_bc) - 1.645*e(se_tau_rb) in `z'
		replace up95_`y'= e(tau_bc) + 1.96*e(se_tau_rb) in `z'
		replace low95_`y' = e(tau_bc) - 1.96*e(se_tau_rb) in `z'
		replace coef_`y' = e(tau_bc) in `z++'
	}

	* Graph
	twoway (rarea up95_`y' low95_`y' r_`y', lcolor(ltbluishgray%8) ///
				msize(*0.7) color(blue%8)) ///
		   (rarea up90_`y' low90_`y' r_`y', lcolor(ltbluishgray%10) ///
				msize(*1.3) color(blue%10)) ///
		   (connected coef_`y' r_`y', mcolor(gs9) lcolor(blue) ///
				msize(medium) msymbol(circle)) ///		
		, yline(0, lpattern(dash) lcolor(red)) ylabel(`yl', nogrid) ///
			graphregion(fcolor(white) color(white)) /// 
			title(`ti', color(black) ) xla(`int') xsize(7) /// 
			legend(off) yti(Estimated discontinuity) xtitle(`xti') ///
			name(r_`y'`gg', replace) nodraw

	drop up90_* low90_* up95_* low95_* coef_* r_*
}


* Final figure
graph combine r_`y' r_`y'l r_`y'h, r(1) graphregion(color(white)) ///
	xsize(12) iscale(1) name(`y', replace)	


*******************************************************************************
* End of Do file
*******************************************************************************
